# Loading in Packages
library(tidyverse)
library(ggtext)
here::i_am("notebook/EEG_ROPE_Sims.Rmd")
here::here()
## [1] "D:/source/repo/psy-project"
# loading in User-Defined Functions
source(here::here("my-scripts/R", "eeg_rope_pipeline.R"))
source(here::here("my-scripts/R", "eeg_sim_functions.R"))
In this step, we use 1/F noise generated from a functions derived from Hank Steven’s work (original function is available on GitHub). After define to templates, one containing an effect and one without, we combine these with noise in order to represent a single trial. We then regenerate the noise again for each participant.
set.seed(314)
.true_onset <- 450 ## Effect at 450ms
.stim_on <- 300 ## Stimulus is shown at 300ms
.freq <- 4 ## Defines resolution of stimulation, frequency of electrode (This is around 250 reading a second)
.max_time <- 800 ## End time for Epoch
.Xf <- seq(0, .max_time, .freq) ## start from zero and end at max time, going up in 4ms (freq)
Nf <- length(.Xf) ## Number of timepoints
# Template one only contains zeros
.temp1 <- vector(mode = "numeric", length = Nf)
# define magnitude of large effect
.magnit <- 2
# Generate ERP Peak
.erp <- dnorm(seq(-1.5, 1.5, length.out= 200/.freq), 0, 1)
.erp <- .erp - min(.erp)
.erp_small <- .erp / max(.erp)
.erp_large <- .magnit*.erp_small
# Getting length and storing them for later use
.l_erp <- length(.erp)
.l_pre_stim <- ceiling(.true_onset / .freq)
# Template 2 contains the ERP peak
temp1 <- c(rep(0, .l_pre_stim), .erp_small, rep(0, (Nf - .l_erp - .l_pre_stim)))
temp2 <- c(rep(0, .l_pre_stim), .erp_large, rep(0, (Nf - .l_erp - .l_pre_stim)))
tibble(x = .Xf, y = temp2) |>
ggplot(aes(x, y)) +
geom_line()+
geom_line(aes(y = temp1))+
annotate("rect", xmin = .true_onset, xmax = (.true_onset+200)-.freq, ymin = 0, ymax = .magnit,
alpha = .5, fill = "forestgreen")+
annotate("rect", xmax = .true_onset, xmin = 0, ymin = 0, ymax = .magnit, alpha = .3)+
annotate("rect", xmax = .max_time, xmin = (.true_onset+200)-.freq, ymin = 0, ymax = .magnit, alpha = .3)+
annotate("text", label = "Null Hypothesis", x = 10, y = .magnit * .97, hjust = 0)+
annotate("text", label = "Alternative Hypothesis", x = .true_onset, y = .magnit * .97, hjust = 0, colour = "forestgreen")+
theme_minimal()+
scale_y_continuous("Time (ms)", expand = c(0,0))
sample_size <-c(10, 25, 50, 100, 150, 200) # number of trials
.gsp <- 1 # gamma spectral power
.outvar <- 1 # noise variance
alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE
effect_time <- c(.true_onset - .stim_on, .true_onset + 200 - .freq - .stim_on)
time_cond <- expression((time < effect_time[1] | time > effect_time[2]))
## storing values in list, so make functions neater
eeg_pipeline_attr <- list(
sample_size = sample_size, # Defines Sample Size/ Number of Trials
alpha = alpha,
num_time_points = Nf,
max_time = .max_time,
cond1_base = temp1,
cond2_base = temp2,
gamma_spec_power = .gsp,
noise_var = .outvar, # Influence Noise output - Zero means no Noise
stim_onset = .stim_on # Time when Stimulus was "shown"
)
## No effect Present (Only using temp1)
eeg_pipeline_attr2 <- eeg_pipeline_attr
eeg_pipeline_attr2$cond2_base <- eeg_pipeline_attr$cond1_base
## Effect Present in template 2 (temp2)
sim_df_1 <- sim_eeg_pipeline(eeg_pipeline_attr, seed = seed)
sim_df_2 <- sim_eeg_pipeline(eeg_pipeline_attr2, seed = seed)
## df contains p, t, ci and rope min and max values, will also be used to plot line of significance at bottom of graphs
rope_df_1 <- sim_df_1 |>
generate_rope_df(eeg_pipeline_attr,
rope_func = get_eeg_nu_rope,
static_margins = static_margins,
replication = F,
sequential = F,
null_rope = T)
##
## A Null ROPE was used to Define Significance
rope_df_2 <- sim_df_2 |>
generate_rope_df(eeg_pipeline_attr2,
rope_func = get_eeg_nu_rope,
static_margins = static_margins,
replication = F,
sequential = F,
null_rope = T)
##
## A Null ROPE was used to Define Significance
.p1 <- sim_df_1 |>
plot_sim_base(rope_df = rope_df_1)
.p2 <- sim_df_2 |>
plot_sim_base(rope_df = rope_df_2)
.p1
## [[1]]
.p2
## [[1]]
# ggsave(plot = p1, filename = here::here("images", "ch4_rope_eff_plain.png"), width = 10, height = 8, dpi = 360)
# ggsave(plot = p1, filename = here::here("images", "ch4_rope_eff_plain.png"), width = 10, height = 8, dpi = 360)
We are going to define two ROPE (Regions of Practical Equivalence) using the data to represent the null and alternative hypothesis. The Null ROPE will be defined by the 90% quantile range of the mean pre-stimulus electrical activity for both groups and the Alternative ROPE will be defined by the mean differences of only the significant p-values. It may be better specify the range of time points a priori, however, I will use two different ROPEs to represent the alternative hypthesis, one defined from raw p-values and the other by bonferroni corrected p-values. This is provide a more conservative definition of an effect, as some results pre-stimulus can be significant and those false positive should not influence the alternative ROPE.
Once the ROPE is created for an “experiment”, it will superimpose the data of the next experiment, to see how the new data fits to the previous study. If the Confidence intervals are perfectly contained in a ROPE, then we can consider that as support for rejecting and maintaining the null hypothesis. For example, if the Confidence interval is inside the Null ROPE, even if the difference is statistically significant, this should likely be interpreted as evidence for keeping the null. Likewise, if the alternative ROPE covers the majority of the ERP amplitude, then that would be evidence of a replication of the underlying effect and thus, evidence to reject the null.
.ch4_sim_1 <- plot_rope(.p1, rope_df_1, null_rope = T, change_subtitle = F)
.ch4_sim_1
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
As sample size increases the effect becomes more clear and obvious, and
the background noise appears to cancel itself out. This in turn
increases the size of the alternative ROPEs, both ROPE defined from raw
p-values and bonferroni corrected are wided at the larger sample size.
However, raw p-values cause the ROPE to cover the entire range of
values, including values where an effect is unlikely. Moreover, it
appears only the Bonferroni defined ROPE maintains a gap between the two
ROPEs, the Null and Alternative, whch would allow easier interpretation.
For, if the ROPEs overlap, and the Confidence interval stretches across
both, then it is unclear which ROPE it supports, preventing us to make a
clear preference of a hypothesis.
It is worth noting that the peak of each of the effects appears to go above the ROPE, as the quantile range which is used to define these ROPEs would need to contain values larger than what is observed. It may not be suitable to use a ROPE to investigate the maximum values, instead it should likely be used to prevent false positives and ensure the observed difference is outside an interval which represents background noise.
To see how a ROPE can prevent false positives, we ran the EEG simulation but with no underlying effect, therefore any statistically significant results are false positives.
.ch4_sim_2 <- plot_rope(.p2, rope_df_2, null_rope = T)
.ch4_sim_2
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
rep_df_1 <- sim_repli_eeg_pipeline(pipeline_attr = eeg_pipeline_attr, seed = seed, total_experiments = 2)
rep_df_2 <- sim_repli_eeg_pipeline(pipeline_attr = eeg_pipeline_attr2, seed = seed, total_experiments = 2)
rep_rope_1 <- rep_df_1 |>
generate_rope_df(pipeline_attr = eeg_pipeline_attr,
rope_func = get_eeg_nu_rope,
static_margins = static_margins,
replication = T,
sequential = F,
null_rope = T)
##
## A Null ROPE was used to Define Significance
rep_rope_2 <- rep_df_2 |>
generate_rope_df(pipeline_attr = eeg_pipeline_attr2,
rope_func = get_eeg_nu_rope,
static_margins = static_margins,
replication = T,
sequential = F,
null_rope = T)
##
## A Null ROPE was used to Define Significance
# ropes <- unique(rep_rope_1$rope)
.rep_p1 <- rep_df_1 |>
plot_rep_base(rep_rope_1)
.rep_p2 <- rep_df_2 |>
plot_rep_base(rep_rope_2)
.rep_p1
## [[1]]
.rep_p2
## [[1]]
.ch4_rep_1 <- plot_rope(.rep_p1, rep_rope_1, null_rope = T, replication = T)
.ch4_rep_1[[1]]
.ch4_rep_2 <- plot_rope(.rep_p2, rep_rope_2, null_rope = T, replication = T)
.ch4_rep_2
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
### Alternative ROPE
rope_df_1a <- sim_df_1 |>
generate_rope_df(eeg_pipeline_attr,
rope_func = get_eeg_alter_rope,
replication = F,
sequential = F,
null_rope = F,
replicate.effect = T)
##
## An Alternative ROPE was used to Define Significance
rope_df_2a <- sim_df_2 |>
generate_rope_df(eeg_pipeline_attr2,
rope_func = get_eeg_alter_rope,
replication = F,
sequential = F,
null_rope = F,
replicate.effect = T)
##
## An Alternative ROPE was used to Define Significance
## adding alternative rope to previous simulations
.ch4_sim_1a <- plot_rope(.ch4_sim_1, df = rope_df_1a, null_rope = F, change_subtitle = T, replication = F)
.ch4_sim_1a
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
## adding alternative rope to previous simulations
.ch4_sim_2a <- plot_rope(.ch4_sim_2, df = rope_df_2a, null_rope = F, change_subtitle = T)
.ch4_sim_2a
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
rep_rope_1a <- rep_df_1 |>
generate_rope_df(pipeline_attr = eeg_pipeline_attr,
rope_func = get_eeg_alter_rope,
replication = T,
sequential = F,
null_rope = F,
replicate.effect = T)
##
## An Alternative ROPE was used to Define Significance
rep_rope_2a <- rep_df_2 |>
generate_rope_df(pipeline_attr = eeg_pipeline_attr2,
rope_func = get_eeg_alter_rope,
replication = T,
sequential = F,
null_rope = F,
replicate.effect = T)
##
## An Alternative ROPE was used to Define Significance
.ch4_rep_1a <- plot_rope(.ch4_rep_1, df = rep_rope_1a, null_rope = F, change_subtitle = T, replication = T)
.ch4_rep_1a
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
.ch4_rep_2a <- plot_rope(.ch4_rep_2, df = rep_rope_2a, null_rope = F, change_subtitle = T, replication = T)
.ch4_rep_2a
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
library(pwr)
df_dyn <- sim_df_1 |>
filter(!eval(time_cond)) |>
group_by(n_trial) |>
summarise(
mean_peak = mean(mn_diff_13),
d = abs(mean_peak)/ sqrt(0.5 * sd (sd_cond1) + sd(sd_cond3)),
half_d = 0.5 * d,
.groups = "drop"
) |>
rowwise() |>
mutate(
beta_2 = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = half_d)$power),
new_alpha_2 = min(0.05, beta_2),
beta = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = d)$power),
new_alpha = min(0.05, beta)
)
dyn_stats <- sim_df_1 |>
pivot_for_stats(replication = F) |>
inner_join(df_dyn, by = "n_trial") |>
mutate(
sig_rule1 = if_else(p < new_alpha, 0, NULL),
sig_rule2 = if_else(p < new_alpha_2, 0, NULL)
)
dyn_stats
## # A tibble: 2,412 × 17
## time n_trial compari…¹ p t diff ci_l ci_u mean_…² d half_d
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -300 10 12 0.287 1.10 0.524 -0.481 1.53 0.374 0.617 0.308
## 2 -300 10 13 0.0419 2.20 0.895 0.0369 1.75 0.374 0.617 0.308
## 3 -296 10 12 0.230 1.24 0.596 -0.415 1.61 0.374 0.617 0.308
## 4 -296 10 13 0.0221 2.50 1.00 0.162 1.85 0.374 0.617 0.308
## 5 -292 10 12 0.228 1.25 0.547 -0.377 1.47 0.374 0.617 0.308
## 6 -292 10 13 0.0345 2.29 0.843 0.0684 1.62 0.374 0.617 0.308
## 7 -288 10 12 0.130 1.59 0.672 -0.219 1.56 0.374 0.617 0.308
## 8 -288 10 13 0.0169 2.64 0.913 0.185 1.64 0.374 0.617 0.308
## 9 -284 10 12 0.184 1.39 0.609 -0.320 1.54 0.374 0.617 0.308
## 10 -284 10 13 0.0246 2.45 0.889 0.128 1.65 0.374 0.617 0.308
## # … with 2,402 more rows, 6 more variables: beta_2 <dbl>, new_alpha_2 <dbl>,
## # beta <dbl>, new_alpha <dbl>, sig_rule1 <dbl>, sig_rule2 <dbl>, and
## # abbreviated variable names ¹comparison, ²mean_peak
.p1[[1]] +
geom_point(data = dyn_stats, aes(x = time, y = sig_rule1 - .9), colour = "darkgreen", inherit.aes = F, size = 1, shape = 15)+
geom_point(data = dyn_stats, aes(x = time, y = sig_rule2 - .8), colour = "darkred", inherit.aes = F, size = 1, shape = 15)
## Warning: Removed 12714 rows containing missing values (`geom_point()`).
## Warning: Removed 2319 rows containing missing values (`geom_point()`).
## Warning: Removed 2275 rows containing missing values (`geom_point()`).
dyn_rep_1 <- rep_df_1 |>
filter(!eval(time_cond), experiment == "A") |>
group_by(n_trial) |>
summarise(
mean_peak = mean(mn_diff_13),
d = abs(mean_peak)/ sqrt(0.5 * sd (sd_cond1) + sd(sd_cond3)),
half_d = 0.5 * d,
.groups = "drop"
) |>
rowwise() |>
mutate(
beta = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = d)$power),
beta_2 = 1 - (pwr.t.test(n = n_trial, sig.level = 0.05, d = half_d)$power),
new_alpha = min(0.05, beta),
new_alpha_2 = min(0.05, beta_2)
)
dyn_rep_stats <- rep_df_1 |>
pivot_for_stats(replication = T) |>
inner_join(dyn_rep_1, by = "n_trial") |>
mutate(
sig_rule1 = if_else(p < new_alpha, 0, NULL),
sig_rule2 = if_else(p < new_alpha_2, 0, NULL)
)
## min difference required to be significant with new alpha rule
min_diff <- convert_t_to_raw(sim_df_1, alpha = dyn_rep_1$new_alpha,
t_values = qt(dyn_rep_1$new_alpha, 100, lower.tail = F), null_rope = F) |>
unlist()
min_diff2 <- convert_t_to_raw(sim_df_1, alpha = dyn_rep_1$new_alpha_2,
t_values = qt(dyn_rep_1$new_alpha_2, 100, lower.tail = F), null_rope = F) |>
unlist()
min_diff_df <- tibble(
n_trial = unique(sim_df_1$n_trial),
min_diff_max = min_diff,
min_diff2_max = min_diff2,
min_diff_min = -min_diff,
min_diff2_min = -min_diff2,
) |>
mutate(across(where(is.double), ~ if_else(!(abs(.) == Inf), ., NULL)))
dyn_rep_stats <- inner_join(dyn_rep_stats, min_diff_df, by = "n_trial")
.rep_p1[[1]] +
geom_point(data = dyn_rep_stats,
aes(x = time, y = sig_rule1 - 1.25, fill = comparison),
inherit.aes = F, size = 1, shape = 22, stroke = NA )+
geom_point(data = dyn_rep_stats, aes(x = time, y = sig_rule2 - 1.5, fill = comparison),
inherit.aes = F, size = 1, shape = 22, stroke = NA )+
geom_ribbon(data = dyn_rep_stats, aes(x = time, ymin = min_diff_min, ymax = min_diff_max), inherit.aes = F, alpha = fill_alpha)+
geom_ribbon(data = dyn_rep_stats, aes(x = time, ymin = min_diff2_min, ymax = min_diff2_max), inherit.aes = F, alpha = fill_alpha)
### Bayes Factors (Testing/ Not Final)
sim_bf_1 <- sim_eeg_pipeline(eeg_pipeline_attr, seed = seed, get_bf = T) |>
mutate(
bf_sig = if_else(bf_13 > 3, 0, NULL),
bf_null_sig = if_else((1/bf_13) > 3, 0, NULL),
)
.p_bf1 <- sim_bf_1 |>
plot_sim_base()
## Argument 'rope_df' was empty - No lines of Significance were plotted
.p_bf1[[1]] +
geom_point(data = sim_bf_1, aes(x = time, y = bf_sig -1 ),shape = 22, fill = fill_colours[1], stroke = NA, inherit.aes = F)+
geom_point(data = sim_bf_1, aes(x = time, y = bf_null_sig - 1),shape = 22, fill = fill_colours[2], stroke = NA, inherit.aes = F)
## Warning: Removed 1006 rows containing missing values (`geom_point()`).
## Warning: Removed 617 rows containing missing values (`geom_point()`).
### Bayesian Parametre Estimation and HDI (Testing/Not Final)
sim_df_1
## # A tibble: 1,206 × 20
## time n_trial mn_co…¹ mn_co…² mn_co…³ sd_co…⁴ sd_co…⁵ sd_co…⁶ mn_di…⁷ mn_di…⁸
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -300 10 -0.184 0.340 0.710 1.00 1.13 0.807 0.524 0.895
## 2 -296 10 -0.326 0.270 0.679 0.920 1.20 0.874 0.596 1.00
## 3 -292 10 -0.278 0.269 0.564 0.815 1.12 0.833 0.547 0.843
## 4 -288 10 -0.304 0.368 0.609 0.848 1.03 0.687 0.672 0.913
## 5 -284 10 -0.290 0.319 0.599 0.835 1.11 0.785 0.609 0.889
## 6 -280 10 -0.233 0.265 0.614 0.755 1.09 0.833 0.498 0.847
## 7 -276 10 -0.328 0.219 0.639 0.815 1.20 0.892 0.546 0.967
## 8 -272 10 -0.441 0.190 0.620 0.850 1.20 0.864 0.631 1.06
## 9 -268 10 -0.538 0.201 0.574 0.681 1.12 0.851 0.739 1.11
## 10 -264 10 -0.527 0.204 0.649 0.640 0.962 0.765 0.731 1.18
## # … with 1,196 more rows, 10 more variables: p_12 <dbl>, p_13 <dbl>,
## # t_12 <dbl>, t_13 <dbl>, ci_l_12 <dbl>, ci_u_12 <dbl>, ci_l_13 <dbl>,
## # ci_u_13 <dbl>, bon_13_sig <dbl>, raw_13_sig <dbl>, and abbreviated variable
## # names ¹mn_cond1, ²mn_cond2, ³mn_cond3, ⁴sd_cond1, ⁵sd_cond2, ⁶sd_cond3,
## # ⁷mn_diff_12, ⁸mn_diff_13
Just like the ROPEs for when there is an underlying effect, the ROPE constructed from raw p-values appears to be less precise and helpful. Here we can see how False positves influence the ROPE whilst no results were significant after bonferroni correction.
The use of a ROPE should reduce the number of false positives, as if the Confidence interval of a statistically significant point is not completely outside the Null ROPE then no conclusions can be made.
It is probably worth replicating this pipeline (the above process) several times to assess behaviour in the long run - as consistency is difficult to assess from a few trials. These following graphs are influenced by randomness, as the noise is randomly generated, so repeating this pipeline should reduce this unwanted influence.
seeds = 1:1000 # the seeds that were used in the 1000 repetitions of the Pipeline
## Reading the Pipeline Simulations - 1000 repetitions
sim_rope_1 <- read_csv(here::here("sim_data", "sim_egg_pipeline_eff.csv"))
sim_rope_2 <- read_csv(here::here("sim_data", "sim_egg_pipeline_no_eff.csv"))
## Evaluating the ROPEs of Number of False Positives
eval_rope_1 <- sim_rope_1 |>
eval_rope(alpha = 0.05, effect_time_vector = effect_time)
eval_rope_2 <- sim_rope_2 |>
eval_rope(alpha = 0.05, effect_time_vector = effect_time)
## Functions to create graphs
source(here::here("my-scripts", "R", "eeg_rope_graphs.R"))
.eval_rope_g1 <- eval_rope_1 |>
fp_graphs(effect = T)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
.eval_rope_g2 <- eval_rope_2 |>
fp_graphs(effect = F)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
Raw p-values are going to generate false positives, which is expected after doing around 400 t-tests. The ROPE does seem to catch some of these, as if the Confidence interval overlaps with the ROPE, then we don’t conclude it is a significant difference and without a decision. However, bonferroni corrected p-values seems to be prefect as controlling for false positives, which is probably because of the multivariate approach used in the analysis. IF the EEG was analysed with the mean amplitude then bonferroni correction will start to resemble the raw p values, and hopefully the ROPE will still be effect.
Whilst minimising False positives is desirable, we can achieve a false positive rate of zero by sacrificing our ability to detect meaningful statistical differences. This is the classic issue of Sensitivity vs Specificity. For example, if we set the margins of a ROPE to -10 and +10, we would have no False Positive because every value is inside that ROPE, so we would have no True positives.
.eval_rope_g3 <- eval_rope_1 |>
tp_graphs()
## [[1]]
## Warning: Removed 3 rows containing missing values (`geom_point()`).
##
## [[2]]
## Warning: Removed 3 rows containing missing values (`geom_point()`).
##
## [[3]]
##
## [[4]]
# auc_df <- eval_rope_df |>
# inner_join(tp_df_1, by = c("n_trial", "rope"))
## Not as neat as I expected - maybe we need more values or the data is not properly formatted?
# eval_rope_1|>
# filter(n_trial != 10) |>
# ggplot(aes(x = fpr, y = tpr, colour = rope))+
# geom_line()+
# theme_minimal()
Next steps: Repeat pipeline to get smooth curves and long-term behaviour of the various methods Repeat pipeline but analyse data using mean amplitude instead of t-tests at every time point Add references to Laken’s and Kruschke’s work, tidy up unprofessional language.
alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE
## defining region where peak occurs - used to find mean amplitude
time_cond <- expression((time < effect_time[1] | time > effect_time[2]))
## Process to proform anovas on each set of n_trial for both conditions
peak_anova_df_1 <- sim_df_1 |>
select(1:4) |> # removing columns that are used for Multivarient Approach
mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |>
pivot_longer(cols = starts_with("mn_cond"), names_to = "cond", values_to = "value") |>
mutate(cond = as.factor(cond)) |>
filter(!eval(time_cond)) |>
nest(data = -n_trial) |>
mutate(model = map(.x = data, ~ anova(lm(value ~ cond, .))),
tidy = map(model, broom::tidy)) |>
select(n_trial, tidy) |>
unnest(cols = c(tidy))
#
# peak_df_1 <- sim_df_1 |>
# select(1:5) |> # removing columns that are used for Multivarient Approach
# mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |>
# group_by(n_trial, in_effect) |>
# summarise(
# mn_peak = mean(mn_diff, na.rm = T),
# mn_peak_t1 = mean(mn_diff, na.rm = T, trim = .2),
# mn_peak_t1 = mean(mn_diff, na.rm = T, trim = .1),
# mn_peak_t5 = mean(mn_diff, na.rm = T, trim = .05),
# sd_peak = sd(mn_diff, na.rm = T),
# .groups = "drop"
# ) |>
# pivot_longer(cols = contains("mn_"),
# names_to = "mean_type",
# values_to = "value")
peak_df_1 <- sim_df_1 |>
select(1:4) |>
mutate(in_effect = if_else(!eval(time_cond), 1, 0)) |>
dplyr::group_by(n_trial, in_effect) |>
dplyr::summarise(
mn_peak_cond1 = mean(mn_cond1),
mn_peak_cond2 = mean(mn_cond2),
mn_diff = mn_peak_cond2 - mn_peak_cond1,
p = t.test(mn_cond1, mn_cond2)[["p.value"]],
ci_l = t.test(mn_cond2, mn_cond1, conf.level = 1 - alpha, alternative = "two.sided")[["conf.int"]][1],
ci_u = t.test(mn_cond2, mn_cond1, conf.level = 1 - alpha, alternative = "two.sided")[["conf.int"]][2],
.groups = "drop"
) |>
dplyr::group_by(n_trial) |>
dplyr::mutate(raw_sig = dplyr::if_else(p < alpha, min(ci_l), NULL)) |> # just so the graph has a line at zero to show significance
dplyr::ungroup()
#
# sim_df_2 |>
# select(1:5) # removing columns that are used for Multivarient Approach
#
peak_df_1
## # A tibble: 12 × 9
## n_trial in_eff…¹ mn_pea…² mn_pea…³ mn_diff p ci_l ci_u raw_sig
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 10 0 -0.0991 0.0409 0.140 3.82e-4 0.0636 0.216 -0.594
## 2 10 1 0.899 0.465 -0.434 5.53e-7 -0.594 -0.274 -0.594
## 3 25 0 -0.0423 0.0156 0.0579 2.91e-3 0.0200 0.0958 -0.317
## 4 25 1 0.723 0.543 -0.180 1.09e-2 -0.317 -0.0423 -0.317
## 5 50 0 -0.00940 -0.0284 -0.0190 8.42e-2 -0.0406 0.00258 NA
## 6 50 1 0.621 0.680 0.0590 3.73e-1 -0.0719 0.190 NA
## 7 100 0 0.00175 0.0115 0.00974 2.99e-1 -0.00870 0.0282 NA
## 8 100 1 0.586 0.556 -0.0302 6.40e-1 -0.158 0.0976 NA
## 9 150 0 0.00403 0.0214 0.0174 2.48e-2 0.00221 0.0325 -0.183
## 10 150 1 0.579 0.525 -0.0538 4.09e-1 -0.183 0.0750 NA
## 11 200 0 0.00558 0.00680 0.00122 8.48e-1 -0.0112 0.0137 NA
## 12 200 1 0.574 0.570 -0.00377 9.55e-1 -0.136 0.128 NA
## # … with abbreviated variable names ¹in_effect, ²mn_peak_cond1, ³mn_peak_cond2
peak_df_1 |>
filter(in_effect == 1) |>
ggplot(aes(n_trial, group = n_trial))+
geom_point(aes(y = mn_peak_cond1), colour = line_colours[1], shape = 15)+
geom_point(aes(y = mn_peak_cond2), colour = line_colours[2], shape = 15)+
geom_hline(aes(yintercept = 0), size = 1.05)+
ggtitle("Mean Amplitude of ERP Peaks in two Conditions")+
labs(subtitle = glue::glue("Conditions:,
<b><span style='color:{line_colours[1]};'>Control Group </span></b>and
<b><span style='color:{line_colours[2]};'>Alternative</span></b>"))+
scale_y_continuous("Mean Amplitude (mV)")+
scale_x_continuous("Sample Size")
# geom_point(aes(y = mn_diff), colour = "grey", shape = 15)
# geom_errorbar(aes(ymin = mn_peak_cond2 - ci_l, ymax = mn_diff + ci_u), colour = line_colours[2])
# ggsave(filename = here::here("images", "ch4_mean_peak_1.png"), width = 10, height = 8, dpi = 360)
peak_df_1 |>
filter(in_effect == 1) |>
ggplot(aes(n_trial, group = n_trial))+
geom_point(aes(y = mn_diff), colour = line_colours[1], shape = 15)+
geom_errorbar(aes(ymin = mn_diff - abs(ci_l), ymax = mn_diff + abs(ci_u)), colour = line_colours[1], width = 5)+
# geom_point(aes(y = mn_peak_cond2), colour = line_colours[2], shape = 15)+
geom_hline(aes(yintercept = 0), size = 1.05)+
ggtitle("Mean Difference in Amplitude of ERP Peaks in Various Sample Sizes")+
scale_y_continuous("Mean Amplitude (mV)")+
scale_x_continuous("Sample Size")+
labs(subtitle = "Error Bars Represent 95% Confidence Intervals of Mean Difference")
# ggsave(filename = here::here("images", "ch4_mean_peak_2.png"), width = 10, height = 8, dpi = 360)